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Abstract 

We consider the Friedmann- Robertson- Walker cosmologies of theories of gravity that generalise 

fv^ ■ the Einstein-Hilbert action by replacing the Ricci scalar, R, with some function, f{R). The general 

asymptotic behaviour of these cosmologies is found, at both early and late times, and the effects of 

. -. ■ adding higher and lower powers of R to the Einstein-Hilbert action is investigated. The assumption 

^h| that the highest powers of R should dominate the Universe's early history, and that the lowest 

I ■ powers should dominate its future is found to be inaccurate. The behaviour of the general solution 

{3JQ| is complicated, and while it can be the case that single powers of R dominate the dynamics at 

late times, it can be either the higher or lower powers that do so. It is also shown that it is often 

the lowest powers of R that dominate at early times, when approach to a bounce or a Tolman 

^ ■ solution are generic possibilities. Various examples are considered, and both vacuum and perfect 

^S| ' fluid solutions investigated. 

00 ■ 

1 Introduction 
!>■ 

O. 

OO , We study here the dynamics of Friedmann-Robertson- Walker (FRW) universes in f{R) theories of 

gravity. These theories are derived from generalisations of the usual Einstein-Hilbert Lagrangian of 

General Relativity (GR), such that 

X' C = f{R), (1) 

d ' and have been considered extensively in the literature (see e.g. [B [2 El S]). Specification of the 

function f{R) defines the theory, and GR can be seen to be the special case f = R. Such theories 
have drawn considerable interest as they found success in early attempts to create a perturbatively 
re-normalisable quantum field theory of gravity [5] , as well as turning up more recently in the effective 
actions of string theory [6l[7]. In cosmology these theories have been used extensively in attempts to 
explain the late-time accelerating expansion of the Universe [U [9] , cosmological infiation |101 [TTl [T2] 
and the nature of the initial singularity [131 [T^ I15j . For a recent review see |16j . 

In considering generalised f{R) theories of gravity it is often implicitly assumed that at late times 
in the evolution of the Universe it should be the lowest powers of R that dominate the gravitational 
Lagrangian. That is, at late times we should have i? — > 0, and the Universe should behave as if it 
were governed by a gravitational Lagrangian of the form 

£o=lim/(ii), (2) 
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which is often presumed to correspond to the Einstein-Hilbert Lagrangian, although other hmits have 
been considered in attempts to address the apparent late-time acceleration of the Universe [51 |9]. 
Conversely, the introduction of higher powers of R into the gravitational Lagrangian has often been 
assumed to mean that at early times the Universe should behave as if governed by the vacuum dynamics 
of a Lagrangian 

£00 = hm f{R). (3) 

R^oo 

The picture is then one of a universe that starts off at high R, dominated by a Lagrangian of the 
form ([3]), and that subsequently expands until R becomes small and the gravitational dynamics are 
dominated by a Lagrangian of the form ([2]). It is the purpose of this paper to determine the veracity of 
such assumptions. This is achieved by studying the dynamical evolution of FRW universes, governed 
by theories with general f{R)- The asymptotic behaviour of the general solution to the Friedmann 
equations is then investigated, and used to evaluate the extent to which the afore mentioned behaviour 
may be considered generic. 

The general solutions of FRW cosmologies governed by f{R) theories of gravity have been studied 
previously by a number of authors, in a number of different contexts. Much of this work has made use 
of the dynamical systems approach, which has been used to study specific classes of f{R) in isotropic 
cosmologies in [171 HSl ES] , and anisotropic cosmologies in [201 IHl ES] . Exact analytic expressions have 
been found for the general FRW solutions of some f{R) theories in [23], and the dynamical systems 
approach applied to general f{R) has been considered in [24)^1 . For studies of spherically symmetric 
and weak field solutions see [Ml ttZl |26l [27], and references therein. The approach used here is a 
generalisation of the analysis performed in [T7|, where theories of the form / oc i?" were considered. 

We find here that the late-time attractors of FRW cosmologies for general f{R) have various 
different forms, and that the asymptote toward which the general solution is attracted depends upon 
the initial conditions. Some of these solutions correspond to the lowest powers of f{R) dominating at 
late-times, and others to the highest powers. Expanding universes with powers of R lower than R^ 
dominating their dynamics generically appear to exhibit the former behaviour, while universes with 
powers of R greater than R^ dominating appear to generically exhibit the latter. Expanding universes 
with higher powers of R dominating their early evolution therefore appear unlikely to evolve to a state 
where the Einstein-Hilbert term dominates. 

We also find that there usually exist multiple early-time attractors for the general solution. These 
can take on different forms, but generically it appears that they either evolve as a ~ ta^ toward a 
big bang singularity in the past, or that they approach a point of inflexion, where the scale factor 
is constant. This is in good agreement with the analytic general solutions for / oc i?" found in [23j . 
These results do not mean that a period of inflation cannot occur (as indeed it appears to if / ~ i?^), 
but it does mean that the general solution does not generally start off inflating (even if / ~ i?^). It 
also means that the picture of the highest powers of R dominating the earliest stages of the Universe's 
evolution may not be an accurate one. 

We begin in section 2 by giving the FRW field equations for f{R) theories, together with some 
simple power-law particular solutions that will later appear as asymptotes of the general solution. 
In section 3 we use a dynamical systems approach to determine the form of the general solution for 
vacuum cosmologies. The phase space of the general solution is two dimensional, and the location and 
stability of critical points in this space is determined. In section 4 we perform a similar analysis for 
the case of perfect fluid cosmologies. The phase space of the general solution is now three dimensional, 
and the location and stability of critical points is again determined. In section 5 we consider the effect 
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of adding higher and lower powers of R to the Einstein-Hilbert action. Section 6 provides a discussion 
of the resuhs, and the appendix gives some special cases that are of particular interest. 

2 f(R) Cosmology 

2.1 Field Equations 

Replacing the Ricci scalar, R, in the Einstein-Hilbert action by a more general function, f{R), gives 
the Lagrangian density 

C = f{R) + Cm, (4) 

where Cm is the Lagrangian density of matter fields. Variation of the corresponding action, with 
respect to the metric, gives the field equations 

fRRab - -^fgab + /iJ;"^ {QabQcd - Qacgbd) = SvrTab, (5) 

where fu = 6f /5R, and Tab is the energy-momentum tensor of matter fields, defined in the usual way. 
The Lagrangian formulation of the theory guarantees the conservation equations T^-^-f, = 0. 

Substituting the FRW metric into the field equations, ([5]), gives the analogue of the Friedmann 
equations 

-3-//? + lf + 3-fR = Sirp (6) 

a 2 a 

-JrR + 2f + 3(fR + S^Jr] = 8^(4 - 3j)p, (7) 

where a = a{t) is the scale factor, and here we have assumed the matter fields are well described by a 
perfect fiuid, with barotropic equation of state p = ("f — l)p- The conservation equations are then, as 
usual, 

p + 3j-p = 0, (8) 

(Jj 

and the Ricci scalar is 

1 fa cP- k\ ,^. 

^=«- + - + -' (9) 

where k is the constant curvature of homogeneous spatial 3-surfaces. Henceforth, unless explicitly 
stated otherwise, we will be considering only the spatially fiat class of models, where k = 0. Spatially 
curved models will be investigated elsewhere. 

2.2 Particular Solutions 

The field equations ([HD-Q have a number of particular solutions that are of interest. These solutions 
often act as attractors, toward which the general solutions can asymptote. Such behaviour will be 
shown in subsequent sections. 

First let us consider the 'vacuum dominated' solution 

a ~ (t — to) ^~^^o+2Sg [vacuum dominated solution] (10) 



where to is a constant, and Bq is defined as 

/ 
Bq = 1 — -— > constant. (11) 

That Bq ^constant is a requirement on the theory in order for the solution to exist. Equation (|1U|) 
represents a cosmology in which the evolution of the scale factor is dominated by the dynamics of the 
Ricci curvature itself. It is this type of evolution that is often invoked to account for the late-time 
acceleration of the Universe, or its inflation at early times. In the limit that GR is approached, and 
Bq -^ 0, this solution reduces to Minkowski space. 

Another solution of interest is the 'matter dominated' solution 



a ~ (t — to)^^*^"^"' • [matter dominated solution] (12) 

Again, this solution exists if Bq — ^constant (not necessarily the same constant as in (llOp ). This 
solution, if it exists, is dependent on the matter content of the space-time, as can be seen from the 
explicit dependence on 7. As Bq -^ 0, and GR is approached, this solution reduces to the usual 
spatially flat, perfect fluid dominated Friedmann solution. 

The last particular solution of interest is the Tolman solution 

a ~ (t — to)^- [Tolman solution] (13) 

This solution is independent of both the matter content of the Universe, and of the form of f{R)- We 
call it the Tolman solution as it is identical to the radiation dominated Friedmann solution of GR. 
However, this solution does not necessarily require the presence of either radiation (7 = 4/3), or of an 
Einstein-Hilbert term in the action. 

We will now proceed to investigate the form of the general solutions to equations Q-Q, first for 
the vacuum case (p = 0), and then for the perfect fluid case (p / 0). The particular solutions above 
will often appear as the asymptotic form of the general solution, as t ^ ±00 or to =constant. 

3 Vacuum Cosmologies 

3.1 The Dynamical System 

In the case of vacuum cosmologies we have p = 0. The field equations ©-Q can now be transformed 
into a system of first-order differential equations, by defining the time coordinate dr = -^/iRldt, and 
the new variables x = a' /a, y = f'^/ fn and z = In \R\. The equations dn])-® then become 

(14) 

(15) 
(16) 

with the constraint 

(17) 

where primes denote differentiation with respect to r, and Q =sign(i?). These definitions ensure r is 
always real, and monotonically increasing in t. The two new functionals A = A{z) and B = B{z) are 
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defined by ^ = fR/2Rfjifi and B = I — f/Rfji. Specifying f{R) gives A{z) and B{z), and the system 
of equations (J14p - (jl7p is closed. However, even before specifying f{R) it is possible to determine some 
generic features of the system above, and hence of vacuum f{R) FRW cosmologies in general. The 
form of equations (J14p - (jl7p allow us to treat the problem as a dynamical system, in which the general 
solutions for f{R) vacuum FRW cosmologies are given as trajectories in the phase space {x,y,z). Such 
an analysis will allow insight into the behaviour of these models. 

First we note that the surface i? = is an invariant sub-manifold of the {x,y,z) phase space. From 
pT|) it can be seen that R = corresponds to x'^ + xy = 0. Equations (fT^ - (fT7|) then give 

{x^ + xy)' = -y{x^ + xy){l + 2A- l/B), (18) 

so there are no trajectories in (x, y,z) that allow R to change sign along them. This justifies our choice 
of z. 

It is now convenient to perform a transformation from the infinite plane spanned by (x, y), to a 
finite closed space. This can be achieved by the re-definitions x = rcosO and y = rsinO, which give 



-Q[2(l-2i3)cos. + (l + i3)sinW ^ + ^°t2'J^+""^'^^ (19) 



z' = 2AP. -^ -, (20) 

Y 6 cot 6^(1 + cot 6*)' ^ ^ 

where the constraint ()17p has been used, P is the sign of y, and 9 runs from to 2it. It can now be 
seen that, independent of the value of B, there exist stationary points in (x,y) at 

0^ = -7r/2, 02 = vr/2, 63 = -n/i and 6^4 = Sir/i. (21) 

In the limit B ^ Bq =constant or itcxD, if this exists, there are two further stationary points in the 
{x,y) plane at 

The value of Bq, and the existence of these last two points, depends on the form of f{R), and can 
be straight-forwardly deduced once this function is specified. Various cases will be considered in later 
sections. First we will investigate the behaviour of the scale factor at the points 1 — 6, identified above. 

From the definitions of x and y we have, at stationary 9i, that 

7^=tan0i-, (23) 

JR a 

which can be integrated to obtain fn oc a*^"^^*. Eliminating / in equations ([6]) and ([7]), and substituting 
for fji, then gives 

Oi ^ {2 + 2tanei-tan'^ei) af 

ai~ (2 + tan0i) af ' ^ ' 

which, for tan^j 7^ or 1, integrates to 

Oj = ao|t-tol *"""»''""'''"'', (25) 
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Table 1: The evolution of the scale factor, a(t\ at the critical points 1-6 in vacuum f{R) cosmologies. 
The listed symbols correspond to those used in Figure [TJ 

where oq and to are constants. Substituting Qi from points 1 — 6 then gives the form of the scale factor 
as these points are approached. For the special cases tan 0j = or 1 we instead have 



aoe 



c(i-io) 



(26) 



where c is a constant. For points 5 and 6 this corresponds to Bq = 1/2 or 1. 

The stationary points 1 and 2 can now be seen to correspond to a! ja = 0, and points 3 and 4 to the 
Tolman solution, ()13p . For B ^ 1/2 or 1 we find points 5 and 6 correspond to the vacuum dominated 
solution, (jlOp . These results are summarised in Table [H where symbols refer to those used in Figure 

m 



3.2 Stability Properties 

Having determined the location of the stationary points in the {x,y) plane, and the scale factors 
to which they correspond, we will now investigate their stability properties. This will allow us to 
determine which points are stable, and can act as asymptotic attractors of the general solution, and 
which are unstable, and can act as repellors (i.e. attractors in the asymptotic past, if time is run 
backwards). 

The stability of points 1 to 6 can be determined as follows. Begin by perturbing 9 as 



where 69i is small. For points at finite z = Zi also perturb z, A and B as 

z ^ Zi + 5zi 

5zi 



B^B(z.)+{^) 



Szi. 



(27) 

(28) 
(29) 

(30) 



Z = Zi 



For points at infinite z we need only to perturb 6, and to check the sign of z' , in order to determine 
stability. 

Points 1 and 2, to linear order in perturbations, then give 



K. 



{l + B 



1,2 



QB 



1,2 



2Si,2 y 3(1 + cos(20i,2) + sin(20i,2)) 



se,: 



(31) 
(32) 



where the upper sign is for point 1, and the lower sign for point 2. This expression is independent of 
any 5z, and so is vahd for all points at finite or infinite z. 

For points 3 and 4 we have 



'^^^'^ ^ 4^3,4 V 3(l + cos(2e3,4)+sin(2e3,4))^^''' ^^^^ 

~ Sign {±i^^|i-lli} 503,4, (34) 

where the upper sign is for point 3, and the lower for point 4. Again, this is independent of any bz. 
Points 1 to 4 are therefore stable attractors if 56' /59 < 0, and unstable repellors if 56' /59 > 0. 

Now consider the stability of points 5 and 6. In this case a knowledge of the behaviour of z is 
necessary to determine stability in the (x,y) plane. For Bq ^ 1/2 or 1 the Ricci scalar corresponding 
to the vacuum solution (jlOp is given by 

_ Bo{l + Bo){5Bo-l) 

''''''- 6{l-3Bo + 2B^,nt- tor ^ ^ 

"^"n V 1^0(1 + i^o)(5i^o- 1)1 ^"" "°T ^^ 

where the upper branch in the final line is for t > to, and the lower branch for t < tQ. As r — > itoo, 
we then have that R^^q -^ itcxD or 0. These points therefore exist at finite 9 and infinite z. As such 
we can check their stability by perturbing 9, as before, and checking the sign of z' as r ^ ztco. From 
()36p it can be immediately seen that for t > to we have z' < 0, so when 59' /69 < this is a stable 
point in (x,y,z), and corresponds to i? ^ as r — > oo. Similarly, for t > to it can be seen that when 
59' /59 > this is an unstable point, corresponding to i? ^ ±00 as r — > -co. The condition for the 
existence of points 5 and 6, when t > to, can now be seen to be -B ^ -Bo as 2; — > =Fc>o, for stable or 
unstable points, respectively. Similar behaviour can be seen to be true when t < to, with the opposite 
limits taken for z in each case. Substituting 6* ^ ^5 g + 59^ q into ([T9|) now gives 




3(1 + cos(26'o) + sin(2(9o)) 

{i+Bor " 



^""^ 69,,, (37) 



~sign|T^^^|^}<^%6. (38) 

The upper sign here corresponds to 9, G (— 7r/2,7r/2), and the lower sign to 9, € (7r/2, 37r/2). For the 
special cases Bq = 1/2 or 1 the asymptotic value of z is a finite constant. In these cases the stability 
analysis must include a perturbation to z, and the corresponding perturbations this induces in A{z) 
and B{z). We perform this analysis in Appendix lAl For these cases we find that it is possible for 
points 5 or 6 to initially act as attractors, while later acting as saddle points. We call this behaviour 
'semi-stable', and describe it more fully in the appendix. 

The stability properties of the critical points 1 — 6 are summarised in Table [21 Here 'A' stands for 
(stable) Attractor, and 'R' for (unstable) Repellor. 

3.3 The General Solution 

We now know the location and stability properties of all stationary points in the {9,z) phase space. 
The possible trajectories of solutions in this space can therefore be deduced, and are represented in 
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Table 2: The stability properties of the critical points 1-6 in vacuum f{R) cosmologies. 'A' denotes a 
stable Attractor, and 'R' denotes an unstable Repellor. 



the plots of Figure [H which show slices through the {x,y) plane. Trajectories on these diagrams are 
confined to the black circles, and are separated into four sub-spaces by the circles and squares that 
denote 9i^2 and ^3^4, respectively. These four regions are labelled i to iv, and trajectories that begin 
in one region are confined to that region for all r. The triangles correspond to the points 

2(25-1) , , 

tang,= ^j^^^y (39) 

Other than points 1 — 4, these are the only places at which 6' can be momentarily zero. The triangles 
can move, as R evolves, and in the limit B ^ Bq they approach the stationary points 5 and 6. Plot 
[T]^a) corresponds to the cases B < —1 and B > 1/5, when the triangles are in regions i and in. Plot 
Wijo) corresponds to — 1 < i? < 1/5, when they are in regions ii and iv. 

Let us first consider the cases depicted in Figure [I^a): B < —1 and B > 1/5. In this figure the 
arrows show the direction of trajectories. Expanding cosmologies have x > 0, and so are restricted to 
regions i and iv. All trajectories in region iv start at point 3 and end at point 1. The other expanding 
trajectories, in region i, begin at either point 2 or point 3, and all end on point 5, if that point exists. 
If point 5 does not exist then either 9c is in perpetual motion in region i, or it moves out of region i 
into either region ii or iv. In the former case there is then no simple late-time attractor: trajectories 
are constantly evolving toward the moving 9c. In the latter case B enters the interval —1<B< 1/5, 
and the trajectories are simply attracted to either point 3 or point 4, as shown in Figure [T]^b) , and 
described below. The collapsing cosmologies with x < 0, and described by trajectories in regions ii 
and Hi, can be seen to behave in a similar way, but with the direction of time reversed, and with 
points 2 and 3 inter-changed with points 1 and 4. 

Now consider the cases in Figure [TJ^b): —1<B< 1/5. The triangles are now found in regions ii 
and iv. Here the dashed arrows show the direction of trajectories when < B < 1/5, and the dotted 
arrows show them for — 1 < 5 < 0. First consider the former case, with < B < 1/5. Expanding 
cosmologies again have x > 0, and so again must be confined to regions i and iv. Now all trajectories 
in region i begin at point 2 and end at point 3. Trajectories in region iv are attracted to point 1 or 
3 in the future, and in the past can be seen to approach 9c. If the stationary point 5 exists, then 
trajectories in this region will approach it in the past. The collapsing cosmologies in regions ii and Hi, 
can again be seen to behave in a similar fashion to the expanding solutions, but with the direction of 
time reversed. The cases in which —1<B<0, depicted by dotted arrows, have exactly the opposite 
behaviour to the description just given. Past and future attractors must then be interchanged, but 
the description is otherwise the same. 

Let us now consider further the asymptotic form of the scale factor as the stationary points are 
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(a) The {x,y) plane when B < — 1 or i? > 1/5. 
In this case the triangles, corresponding to the 
critical points (|39|l . are located in the intervals 
(— 7r/4,7r/2) and (37r/4,— 7r/2). Arrows denote the 
direction of trajectories around the circle. 
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(b) The (x,y) plane when -1 < B < 1/5. Trian- 
gles now correspond to (|39|) . and are located in 
the intervals (-7r/2,-7r/4) and (7r/2,37r/4). The 
dashed arrows denote the direction of trajectories 
when < B < 1/5, and the dotted ones when 
-1<B<0. 



Figure 1: The {x,y) plane in the phase space of vacuum solutions, defined by equations p^ - (fT6]) . 
The arrows indicate the trajectories of solutions, confined to the circle in each plot by the constraint 
equation, pT|) . Stationary points are indicated by dots, squares and triangles, and divide the phase 
space into the four sub-spaces i to iv, between which trajectories cannot pass. Dots correspond to 
the points 1 and 2, and do not move. Squares correspond to points 3 and 4, and also do not move. 
The Triangles correspond to the critical points defined by (j39p . which can move. They approach the 
stationary points 5 and 6 in the limit B ^ Bq =constant. 



approached. The functional form of the scale-factor, 
above, and is given by (fTOj) . If 

^5(1 + ^5) 



a{t), at points 5 and 6 has been determined 



(1 - W5 + 2Bi] 



>0 



(40) 



then the late-time behaviour of the expanding point 5 approaches a simple scaling solution as i ^ 00. 
If this condition on B^ is not met, then in order to maintain H > it must be the case that t < to- 
In the limit that t ^ to the scale factor then diverges, and a big-rip singularity occurs. It can be 
seen that divergent behaviour occurs if — 1 < B5 < 0, or 1/2 < i?5 < 1. This range of B contains 
some interesting f{R) theories, as we will show later on. Again, the behaviour of a{t) as point 6 is 
approached is similar, but with the direction of time reversed. 

The points 3 and 4 have a simpler interpretation, as, independent of B, the scale factor evolves like 
a radiation dominated flat Friedmann universe when they are approached. Point 3 always corresponds 
to an expanding universe. For < i? < 1/5 this expanding solution is stable, and trajectories that are 
attracted toward it then approach the simple scaling behaviour (J13l) as t ^ 00. For B < oi B > 1/5 
the expanding solution at 3 is a repellor. Trajectories that originate from this point therefore have a 
big-bang singularity in their past, which they approach as a oc (t — to)^ when f — > to- The collapsing 
trajectories corresponding to point 4 have a similar interpretation, with the direction of time reversed. 



Finally, consider the points 1 and 2, where x/y -^ and 1/y -^ 0. These points correspond to 
a' /a — >0, which is either a stationary point in the evolution of a{t), an approach to Minkowski space 
or a divergence in R. In order to determine which of these is the case it appears necessary to have a 
knowledge of the functional form of f{R), so that equations (J14p -()17p can be solved for in these limits. 
We can, however, make some progress in the cases in which Ai^2 and Bi^2 ^constant. Equations (llSp . 
(|16p and (J17p can then be integrated to give 

X ^ ^QB(l + A)(t - to) (41) 

2A 

^^^o + 7^-^p^ln(r-ro), (43) 

where tq and zq are constants. The definition dr = y^|i?|(it can also be integrated to give 

(r-ro)oc(t-to)('+^). (44) 

It can now be seen that y -^ ±00 corresponds to r — > tq. For A > —1 this corresponds to t ^ to, 
while for A < —1 it corresponds to t — > ±00. The scale factor is now found to be 

ln(-')=(r-ro)V(t-to)'(i+^), (45) 

where oq is a constant. A solution that approaches 61^2 with A > —1 reaches a stationary point in 
its expansion at a finite time, to- This solution must then be matched onto another solution with a 
similar stationary point in its past. If j4 < —1 as 81^2 is approached, then these points are reached 
only as t — > 00. These trajectories then approach Minkowski space asymptotically. 

4 Perfect Fluid Cosmologies 

4.1 The Dynamical System 

Now consider the case p ^ in equations dnj-Q. In this case we begin by transforming the time 
coordinate from t to T via 

.T.^*. ,46) 

where it has been assumed p ^ 0, so that T is always real and increases monotonically in t. Defining 
the new variables w = y^\fji\\R\/8-Kp, x = a' /a, y = f'ji/fR and z = ln\R\ the equations ([6|)-([9]) can 
then be written as 

w' = w((^ + Ajy + ^jxj (47) 

x' = Q7 -2\x^ + ^xy + ^Qw'' (48) 

y' = (^7 + l]xy- iy2 _ Q^2 _ 3^^ ^ 4^2 ^49) 

/ = 2Ay, (50) 

10 



with the constraint equation 

x'^ + xy - -QBw'^ - S = 0, (51) 

where primes are now understood to denote differentiation with respect to T, and Q and S are the 
signs of R and Jr, respectively. The functionals A = A{z), and B = B[z) are the same as before. 

It can again be seen that the surfaces i? = are invariant sub-manifolds of the phase space. These 
surfaces are now given by (|5ip as x^ + xy — S = 0, and it can be seen that 

[x^ + xy- Sy = (x^ + xy- S){3'yx + y/B), (52) 

showing that there exist no trajectories in this space along which R changes sign, justifying our choice 
of z. 

To find the stationary points, at finite distances in x and y, we can eliminate w using the constraint 
equation, ([5T]) . The dynamical system is then described by a three-dimensional closed system of 
equations, with four stationary points at finite distance in the {x,y) plane. These points are at 

±^ , (4-37) /^ .„. 

and 

±V2S t3jBoV2S .^. 

xq 10 = — , yq 10 = — , 54 

where the it signs should be chosen consistently for each point, and where B ^> Bq =constant. If the 
limit B ^ Bq =constant does not occur, then points 9 and 10 do not exist. These points must also 
have real Xi and yi in order to exist. 

We will now find the form of the scale factor for the points 7-10. From the definitions of x and y it 
is clear that at x = Xi and y = yi we should have Jr oc e^*^ °^, and a oc e^^^ °K The conservation 
equation ([8|) then yields p oc a^^'^ oc e-37a;i(T-To)^ which, for B ^ 1, allows (j46p to be integrated to 

(t-to)ocei(^'+37x.){T-To)_ (55) 

The scale factor can then be written, as a function of t, in the power law form 



Oj oc (t-to)='^+«''"^- (56) 

At points 7 and 8 the scale factor can be seen to evolve like the Tolman solution, (J13p . while at points 
9 and 10 it evolves like the matter dominated solution, (J12p . One may expect a different behaviour at 
these points when B = 1, but for this value of B points 9 and 10 do not exist for any positive 7. 

Now consider stationary points at infinite distances in the {x,y) plane. To identify these points it 
is convenient to transform to polar coordinates via x = r cos (j) and y = r sin </>. Re-defining the radial 
coordinate as r = r/(l — r) maps the infinite (x,y) plane to a unit disk, where r — > 1 in the limit 
f — >■ 00. In these coordinates equations (j48p and (j49p become 



r 



(67 — 5)r COS0 — 3cos(3(/))r — 3(47^(1 — r) — r )snicj) 

+ 5r2 sin(3(/>) + 2Q(1 - rf {cos (/> - 2 sin (t))w'^ 
cos 0(3r2 - 675(1 - rf + hr^ cos(20) 

+ 3r2 sin(2(/>)) - Q(l - rf(2 cos 4> + sin(/>)u>2 
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2r(l - r) L 



where w is given by (j5ip as 



2 _ (r2 - 25(1 - r)2 + r^ cos(20) + r^ sin(2(/>)) 
"" " SQ(l-r)2 • ^^^> 

As r — > 1, and infinite distances in the {x,y) plane are approached, these equations become 

"' ^(l-(5-67)S)cos0 + 3(l-5)cos(3(/)) -(l-35)sin(/)- (l-5S)sin(30) (58) 

1 + cos(20) + sin(20)) r,^ ^ ^^ ^.^ ^ ^ 2(1 - 25) cos J , (59) 



r 

45 



25(1 -r) 

independent of the behaviour of B. It can now be seen from (159p that stationary points at infinite 
distances, as r ^ 1, can exist only if 1 + cos{2<j)) + sin(20) = 0, or tan(j) = 2(25 — l)/(5 + 1). This 
gives stationary points at 

</.i = -7r/2, (/)2=vr/2, ./.g = -7r/4, (/'4 = 37r/4, (60) 

and at the two solutions of 

tan</>5,6 = 2(25o-l)/(5o + l) (61) 

that exist in the range (0,27r). These points can be seen to be in similar positions to those in the vacuum 
cosmologies, and we will see that the scale factor evolves in a similar way as they are approached. 
Points 1 to 4 exist for any 5, while 5 and 6 require B ^ Bq =constant as they are approached, in 
order to be stationary. 

To find the evolution of the scale factor at points 3-6 we first see that the definitions of x and y 
allow an integral such that //j oc a*'^"''^', just as in the vacuum FRW case. Substituting this into (j48p 
and ()49p allows us to obtain 

a/[ ^ ((37 + 2) + 1(2 + 7) tan 0, - i tan^ 0i) g'^ _ 3-fS 
a {2 + tan(j)i) c? (2-|-tan0j) 

as these points are approached. This equation can then be integrated to give 

X^ = ^r = -^ ■ ^— + Cia ' ^ (2+tan0,) (53) 

a? 67 + (4 + 37) tan (t)i - tan^ (pi ^ ' 

^^^^37+tan0,{|^i^^ (64) 

where the second line is in the limit x^ — > oo, as is the case for all of points 3-6. This equation can be 
integrated again, for tan (pi / or 1, to find the form of the scale factor, as points 3-6 are approached, 
to be 

2(2+tan0j) 

a (X {T - To)~('5^+(*+3'^)'™"*»-*'"'^'^i) (65) 

2cot <^i + l 

OC (t-to) t="i^.-i , (66) 

where Tq and to ^ire constants of integration, and the second line has been obtained by an integration 
of equation (jl6]) . the definition of T. Substituting (pi into the expression above shows that points 3 
and 4 correspond to the Tolman solution, (fT3]) . and points 5 and 6 to the vacuum dominated solution, 
(jlOp . This is as may have been expected by analogy to the vacuum cosmologies, investigated above. 
The special cases taiKpi = and 1, corresponding to points 5 and 6 when Bq = 1/2 or 1, give 
a^fi oc exp{c(t — to)}- Points 1 and 2 correspond to a' /a = 0. 

The evolution of the scale factor at the stationary points 1 to 10 is summarised in Table 3. The 
symbols in this table refer to those used in Figures [2][6l 
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Point 


a(t) 


Symbol 


1,2 


d ja = 


Circle 


3,4 


Tolman solution. (jlSp 


Square 


5, 6 


Vacuum dominated solution. (|10p 


Triangle 


7,8 


Tolman solution. (jlSp 


Square 


9, 10 


Matter dominated solution, (I12p 


Star 



Table 3: The evolution of the scale factor, a(t)^ at the critical points 1-10 in perfect fluid filled /(-R) 
cosmologies. The listed symbols correspond to those used in Figures [2][6j 



4.2 Stability Properties 

Having found all stationary points, at both finite and infinite distances in the (a;,y) plane, we will now 
proceed to establish their stability properties. This will allow some insight into the degree to which 
they can be considered the asymptotic limits of the general solutions to the Friedmann equations, ([6]) 
and ©. 

First consider the points at finite distances in {x,y), points 7-10. As these points are approached 
z can be seen to diverge to ±00. We can therefore check stability in the (x,y) plane by perturbing x 
and y as 

X ^ Xi + u and y ^ yi + v, (67) 

and checking the signs of the eigenvalues, Aj, of the linearised equations 



u' = XiU and v' 

Substituting ([67|) into (|i8]) and ([Ml) gives the linearised system 



XiV. 



{2{2xi + yi) + {{6j-8)xi + y,)B) Xi{2 + B) 

U -\ ^t; V 



2B 

{{IQxj + 2y, + ?,-iyi)B - 4(2x, + y,)) 

2B 



u + 



2B 

(((2 + 37)x. 



2yi)B- Ax, 



2B 



(68) 

(69) 
(70) 



where B in these expressions is in the limit i? — > or 00, whichever is appropriate. The eigenvalues, 
Aj, are then given by the roots of a quadratic equation of the form A? + aXi + /? = 0. If a > and 
/3 > we have a stable point in the {x,y) plane. If a < and /? > we have an unstable point, and if 
/? < we have a saddle point. The behaviour of z can be found by checking the sign of z' . For points 
7 and 8 we find that 



a 



T 



{{4 - 3-f) - (5 - 6j)B)^ 
B^5 - 37 



and 



/3 



(4-37(l-i?))S 



B 



(71) 



where the branch of =F should be chosen consistently with equation ()53p . Similarly, for points 9 and 
10 we find 



a 



3((1 + B)7- 2)^/5 
V2(2-i?(4 + 37(l + i?))) 



and 



P 



(i-3jil-B))S 
B 



(72) 



with ^ chosen consistent with (|54p . It can immediately be seen that one of these two sets of points 
always corresponds to a pair of saddles, as /? has opposite signs in (|7ip and (j72p . The other set can 
then be seen to contain one attractor and one repellor, as a has a different sign for each point in both 
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(|7ip and (j72p . Which two points are saddles, which is the attractor and which is the repellor depends 
on the values of -B, 7 and S. We will consider various cases below. 

Now consider the stability of the stationary points at infinite distances in {x,y). In this case we 
perturb (p so that (p ^ (pi + Sep. For points 1 and 2 we must also consider a perturbation to z such 
that 2; — > Zj + 6z, and B ^ B + {dB/dz)6z. The evolution equation ([58]) and ([59|) then become 



(^3^5</<i,2 and r,,2 = i^' 



1,2 = ^Ft:^ IT^<^</'1,2 and r, 2 = ±7T> (73) 



independent of 5z. The upper sign here is for point 1, and the lower sign for point 2. These points 
are considered stable if r' > and 5(l)/(p < 0, unstable if r' < and 5(p/(p > 0, and as saddle points 
otherwise. Therefore, if i? > or < — 1 then one of these points is stable and the other unstable. If 
— 1 < i? < then both points are saddles. 

For points 3-6 both r and z — > co, while (p is finite (for B ^ 1/2 or 1). It is therefore sufficient to 
perturb (p as (p ^ (pi + 6(p, and to check the signs of r and z. As (^ ^ (pi, and r — > 00, we then have 
for points 03^4 



V2B (1-r) ^^•^ ''''-^^7T 



3,4 = T^^^7r^ and r^_^ = ^^__^, (74) 



where the upper sign is for 03, and the lower sign for (f)4. For points 5 and 6 

A^' ^ (1 - 5B) <^05,6 , , ^ (2 -(4 + 37(1 + i?))i?) ,„. 

^l + TT^^ ^ ^ 2^(1 + ^)^1 + 1^^^ 

where the upper sign is for point 5, and the lower for point 6. It can be seen from the 5(p equations 
that for < B < 1/5 the only points that can be stable are 3 and 6, and that the only points that 
can be unstable are 4 and 5. For B < or B > 1/5 this behaviour is reversed and 4 and 5 are 
the only points that can be stable, while 3 and 6 are the only unstable points. The special cases 
B = 1/2 and 1 can have z asymptoting to a finite value, and so in this case the stability analysis must 
include a perturbation to z. In Appendix [A] we perform this analysis, and find similar behaviour to 
the corresponding vacuum cosmologies investigated above. 

The stability of the critical points 1-10 can be seen to be strongly dependent on the value of 
the parameter B, as in the vacuum case. Now, however, the additional parameters 7 and S are 
also involved in determining stability. The results for the different cases 5 = ±1, and 7 ^ 5/3 are 
summarised in Tables SHTl Having determined the stability properties of all stationary points, at both 
finite and infinite distances in the (x,y) plane, we will now proceed to investigate what this can tell 
us about the behaviour of the general solution. 

4.3 The General Solution 

The general solutions to the Friedmann equations ([6])- ([7]) in the presence of a perfect fiuid are more 
complicated than the vacuum case. The phase space of solutions is now, in general, three dimensional. 
Nevertheless, it is still possible to make progress in understanding the general solution. 

As was previously noted, the surfaces x'^ + xy = S (corresponding to i? = 0) are invariant sub- 
manifolds of the phase space. This equation has two roots, describing two non-intersecting surfaces 
that separate the the phase space into three different regions, between which trajectories cannot cross. 
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Point 


Stability when S 


= — 1 and 7 > 


5/3 


B < -1 


-KB <0 


< B < 1/5 


1/5 <B 


1 


A 


S 


A 


A 


2 


R 


s 


R 


R 


3 


S 


s 


A 


S 


4 


S 


s 


R 


S 


5 


m < 0: A 
■m > 0: S 


A 


m< 0: S 
m > 0: R 


A 


6 


m < 0: R 
■m > 0: S 


R 


m< 0: S 
m > 0: A 


R 


7 


R 


R 


S 


n < 0: S 
n > 0: R 


8 


A 


A 


S 


n < 0: S 
n > 0: A 


9 


m < 0: S 
m > 0: - 


- 


m < 0: R 
m > 0: - 


n < 0: R 
n> 0: S 


10 


m < 0: S 
m > 0: - 


- 


m < 0: A 
m > 0: - 


n. < 0: A 
n> 0: S 



Table 4: The stabihty properties of the critical points 1-10 in perfect fluid f{R) cosmologies, when 
S = —1 and 7 > 5/3. Here 'A' denotes a stable Attractor, 'R' an unstable Repellor and 'S' a Saddle 
point. The quantities m and n are defined as m = 2 — 37-B^ — B{A + 37), and n = 3jB^ + i?(4 — 37). 
Dashes indicate that a point is not present for that range of B. 
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— 1 and < 7 
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-KB <0 


< B < 1/5 
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m > 0: S 


A 
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m < 0: A 
m > 0: S 


6 

7 


m < 0: R 
m > 0: S 


R 


A 


m < 0: R 
m > 0: S 


8 
9 


■m < 0: S 
m > 0: - 





- 


m< 0: S 
m > 0: - 


10 


m < 0: S 
m > 0: - 


- 


- 


m< 0: S 
m > 0: - 



Table 5: The stability properties of the critical points 1-10 in perfect fluid f{R) cosmologies, when 
S = —1 and < 7 < 5/3. Here A, R, S, and m are defined as in Table [H Dashes again indicate the 
absence of a point. 
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Stability when S = 1 and 7 > 5/3 
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-1 
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< B < 1/5 


1/5 <B 
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A 


2 
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4 
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A 


m< 0: S 
■m > 0: R 


A 


6 
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m < 0: 
m > 0: 


R 

S 


R 


m< 0: S 
m > 0: A 


R 


8 
9 


m < 0: 
m > 0: 


A 


At 


m < 0: - 
m> 0: S 


— 


10 


m < 0: 
m > 0: 


R 


Rt 


m < 0: - 
m > 0: S 


- 



Table 6: The stability properties of the critical points 1-10 in perfect fluid f{R) cosmologies, when 
S = 1 and 7 > 5/3. A, R, S, and m are defined as in Tabled! Dashes indicate the absence of a point. 
The t here indicates that these points have the indicated stability as long as 7 < 2. The * indicates 
the range of B for which, in this case, there exist regions of the phase space without either an attractor 
or a repellor. 
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Stability when S = 1 and < 7 < 5/3 
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R 
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A 
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m > 0, ra > 0: S 


ra < 0: A 
ra > 0: S 


ra < 0: A 
ra > 0: S 
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m > 0, ra > 0: S 


ra < 0: R 
ra > 0: S 


ra < 0: R 
ra > 0: S 
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ra <0: S 
ra > 0: A 


ra < 0: S 
ra > 0: A 


ra <0: S 
ra > 0: A 


rai< 0: - 
m > 0: A 
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ra <0: S 
ra > 0: R 


ra < 0: S 
ra> 0: R 


ra <0: S 
ra > 0: R 


m < 0: - 
rai > 0: R 



Table 7: The stability properties of the critical points 1-10 in perfect fluid f{R) cosmologies, when 
S = 1 and < 7 < 5/3. A, R, S, and m are defined as in Table HI Dashes indicate the absence of a 
point. 
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The shape of these regions depends on the sign of S, and is illustrated in Figures [2]l6] by the grey 
and white areas. These plots show the {x,y) plane, and some representative trajectories within it, for 
various values of 7 and B. The left-hand plot, (a), in each of these figures corresponds to 5" = — 1, 
and the right-hand plot, (b), to S = 1. 





(a)S 



(b)S: 



Figure 2: The (x,y) plane in the phase space of solutions defined by equations (|Tr|) - ([5T]) . Here the 
infinite {x,y) plane has been compacted to a unit disk. Dots, squares, triangles and stars indicate 
the position of stationary points. Squares correspond to points 3, 4, 7 and 8, and hence the Tolman 
solution, (fT3l) . Stars correspond to points 9 and 10, and the matter dominated solution, (fT2|) . Triangles 
are at points 5 and 6, corresponding to the vacuum dominated solution, ([10]). Dots are at points 1 
and 2, where a' /a = 0. The phase plane is separated into three different regions by the invariant 
sub-manifolds where R = 0. Grey regions have R < 0, and white regions have R > 0. In these 
plots we take B = 2/3 and 7 = 1, corresponding to pressureless dust and a gravitational Lagrangian 
dominated by a term ~ R^. 



A number of parallels can be drawn between the perfect fluid cosmologies, illustrated in Figures 
[2]l6l and the vacuum cosmologies, illustrated in Figure [H The circle at r = 1 in Figures [2][6] can be 
seen to be quite similar to the circle in Figured! They both have stationary points at tt/2 and —it/2, 
corresponding to a' /a =constant, and at — 7r/4 and 37r/4, corresponding to a ~ ta (although the 
coordinates have different definitions in each case). Whatsmore, in both cases these points mark the 
boundaries of regions between which trajectories cannot pass. A further similarity is in the position of 
the triangles, which correspond to the vacuum dominated solution, (jlOp . In both cases, for B > 1/5 
or < —1, these points exist in what are labelled as region i and in in Figured! and for —1 < B < 1/5, 
they exist in regions ii and iv. 

All of the points at r = 1 therefore have analogous points in the vacuum cosmology case, and all 
of these points can be seen to exist for any 7, B and S (as long as -B — > i?o =constant in the case 
of points 5 and 6). As r increases these perfect fluid cosmologies therefore approach the behaviour 
of the vacuum cosmologies studied above. The principle difference between the perfect fluid case 
and the vacuum case is due to the existence of the extra dimension in the phase space. Points that 
were previously attractors (or repellors) in vacuo, can now be saddles in the presence of a perfect 
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(a) 5: 



{h)S- 



Figure 3: These plots show the {x,y) plane in the phase space of solutions when B = 2 and 7 = 1, 
corresponding to pressureless dust and a gravitational Lagrangian dominated by a term ~ R~^. See 
the caption of Figure [2] for further details. 



fluid, as they can be unstable (or stable) in the extra dimension. Whether these points maintain the 
attractor/repellor nature they exhibited in vacuo, or become saddles, can be read off from Tables |3][71 

As well as the points at r = 1, there can exist four further points at finite distances in the {x^y) 
plane. These are points 7-10, given by equations (f53]) and ([541) . They have no analogy in the vacuum 
cosmologies considered above, and do not exist for all 7, B and S. Points 7 and 8 correspond to the 
Tolman solution, ()13p . and exist on the boundaries between grey and white regions. Points 9 and 10, 
correspond to the matter dominated solution ()12p . and can exist in either the grey or white regions. 
For S = —1 they are saddles if they exist in the the white region, and attractor/repellors if they exist 
in the grey regions. For S = 1 this behaviour is reversed. If both of these sets of points exist then one 
of them will be a pair of saddle points, while the other will be an attractor/repellor pair. In general, 
the existence and stability of these points are functions of 7, B and S, the details of which can be 
read off from Tables HJlTl 

The existence, location and stability of the critical points, in the presence of a perfect fluid, is more 
complicated than the vacuum case. Rather than discuss all possible trajectories, we will therefore show 
a number of examples that we consider to be illustrative. In these examples B will taken as a constant. 
It should be emphasised that in general B is not a constant. However, for the purposes of illustrating 
the existence and attractor/repellor nature of the stationary points, for various different values of B, 
it is convenient to take it as so. 

Figure [2] shows the {x,y) plane for a universe filled with pressureless dust, 7 = 1, for the case 
B = 2/3. Such a value of B may be expected if f{R) ~ B? , and so may be of interest when higher 
powers of the curvature are expected to dominate the gravitational Lagrangian. Here the points 
at r = 1 have locations and stability analogous to the vacuum case. There are also two points at 
finite distance, corresponding to the matter dominated solution, (J12p . when S = —1, and the Tolman 
solution, (113p . when S = 1. These points, however, are both saddles, and asymptotic behaviour in 
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(a) 5: 



{h)S- 



Figure 4: Plots showing the {x,y) plane in the phase space of solutions when B = 1/10 and 7 = 1, 
corresponding to pressureless dust and a gravitational Lagrangian dominated by a term ~ R~9 . See 
the caption of Figure [2] for further details. 



this case is therefore very similar to the corresponding vacuum cosmology. 

Consider now the effect of lower orders of R dominating the gravitational Lagrangian. If f{R) ~ 
1/R dominates then we may expect S — > 2. This case, with a dust equation of state, is illustrated in 
Figure [3l The asymptotes of the trajectories in Figure [3] can be seen to be very similar to those in 
Figure [21 and again the stability and location of the critical points are analogous to the corresponding 
vacuum cosmology. Only saddle points are introduced at finite distances in {x,y). 

Let us now consider some ranges of B for which we know the stability properties of the critical 
points to change. First consider B in the range < i? < 1/5, say B = 1/10. This case would be 
expected from f{R) ~ R~ , and is shown in Figure U for 7 = 1. For S = —1 the points 3 and 4, 
at — 7r/4 and 37r/4, are now no longer attractor/repellors, but saddle points. Trajectories that would 
otherwise have ended on them are now therefore drawn to points 1 and 2, at — 7r/2 and it/2. There 
are no new points at finite distances in this case, and the behaviour is otherwise the same as the 
corresponding vacuum cosmology. For S = 1 points 3 and 4 are again saddles. Now, however, all 
four possible points at finite distances exist. The points 9 and 10 are now attractor/repellors, and 
so this case is qualitatively different from the corresponding vacuum cosmology, with new asymptotes 
possible. The case with 5 = 1 is that which was considered in |17j . 

Another range of B which we may wish to consider is — 1 < -B < 0, which is shown in Figured] 
for B = —1/2, corresponding to f{R) ~ R'^ ■ The stability of the triangular points, corresponding 
to the vacuum solution (jlOp . is now reversed, as would be expected from analogy with the vacuum 
cosmology. Now the circular points at — 7r/2 and 7r/2 are saddles, and the squares again have the 
attractor/repellor behaviour that would be expected from the corresponding vacuum cosmologies. 
The case S = —1 contains no extra critical points at finite distances, and so again has asymptotics 
that are otherwise the same as the vacuum case. For 5" = 1 all four points at finite distances are 
present, and again points 9 and 10 act as attractor/repellors. 
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(a) 5: 



{h)S- 



Figure 5: Plots showing the {x,y) plane in the phase space of solutions when B = —1/2 and 7 = 1, 

2 

corresponding to pressureless dust and a gravitational Lagrangian dominated by a term ~ i?3. See 
the caption of Figure [2] for further details. 



These examples have been intended to show some representative illustrations of the {x,y) phase 
plane, for different values of B. The figures shown above do not show all possible behaviours, and 
should not be considered exhaustive. It is, for example, quite possible to have points 7 and 8 acting 
as attractor/repellors, although this was not explicitly the case in any of the examples above. The 
reader will be able to find ranges of B and 7 where this behaviour occurs, for either S, by referencing 
Tables mil above. 

It is of interest to note that for all but one configurations of i?, 7 and S, there exists in each region 
of the {x,y) plane at least one attractor and one repellor. This is a necessary condition for there to 
exist no closed orbits, and for there to be at least one available attractor for each trajectory to end one, 
and one repellor from which it can begin. The exceptional case is 5 = 1, — l<i?<0 and 7 > 5/3. In 
this case there exist regions which are either without an attractor, or without a repellor. An example 
is shown in Figure [6j The star in the right-hand white region can be seen to act as an attractor, but 
following trajectories in this region backwards shows that they continue to spiral outwards forever. 
Such a trajectories therefore has no simple asymptote in the past. Similarly, in the left-hand white 
region of this plot all trajectories begin on the star, and spiral outwards forever as they are followed 
forward in time. 

Let us now consider further the behaviour of the scale factor, a(t), as the critical points are 
approached. Trajectories approaching the points 1-6 all have evolutions similar to those approaching 
the analogous points in the vacuum case. If condition (j40p is met, then expanding trajectories with 
these points in their past all start with a big-bang singularity. Expanding trajectories that asymptote 
toward these points exhibit power-law behaviour, as t ^ 00. If condition (j40p is not met, then 
expanding trajectories approaching points 5 and 6 have a big-rip singularity in their future. 

The matter dominated points, 9 and 10, have no analogy in the vacuum cosmologies. Trajectories 
ending on these points approach the matter dominated solution (J12p . When B < 1 these points 
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Figure 6: Plots showing the (x,y) plane in the phase space of solutions when B = —1/2, S = 1 and 
7 = 2. This situation corresponds to a stiff fluid and a gravitational Lagrangian dominated by a term 
~ Rs. The left-hand white region in this plot has only a single repellor, and no attractor. Trajectories 
in this region can then be seen to originate from the repellor, but have no stable critical point to attract 
them in the future. They therefore spiral outwards forever, away from the repellor. The right-hand 
white region can similarly be seen to have an attractor but no repellor. See the caption of Figure [2] 
for further details. 



correspond to a big bang in the past of expanding solutions, and to power-law expansion in their 
future. If i? > 1, for expanding solutions, then there is a big rip in the future, where a ^ oo at finite 
t. 

The two remaining points, 1 and 2, are at y ^ itcxD, and x/y -^ 0. In general it is required to 
solve the dynamical equation (J47p - (j5ip . for some particular f{R), in order to find the form of a(t) as 
these points are approached. If A and B ^constant in this limit then (|47p - ()5ip give 



w -^ wo{T - To) 



X 



1+2A 



^S{T-To) + \QBwliT-To)'- 



{T - To) 
zo+4Aln{T-To), 



(76) 

(77) 

(78) 
(79) 



where Tq, wq and zq are constants. The definition of T, ([^6]) . can then be integrated to 

(t-to) = ^(^-To)^ 

2e^ 



(80) 



where to is another constant. The limit y — > oo can now be seen to be reached as T — > Tq, and hence 
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t — > to- The scale factor can then be found, by integrating x, to be 

For B > or B < —1 the scale factor therefore reaches a finite value, aQ, as t — > to- The trajectory 
must then be matched onto another solution at this point, in order to make a complete history. If 
— 1 < -B < then a ^ oo as these points are approached, and there is a big rip as t — > to- 

5 The Effect of Higher Powers 

Most of the asymptotic behaviours found above correspond to scale factors obeying power laws of the 
form a ~ (t — to) , where k is some constant. In such cases the Ricci scalar is given by 

_ 6k{2k - 1) 

Solutions approaching these asymptotes, as t — > ±00 or to, therefore have i? ^ or itcxD. It is then 
straightforward to read off the values of A and B for different f{R)- Consider, for example, a theory 
of the type 

fiR) = Y,CiR\ (84) 

i 

where the q are constants, and the sum is over finite i. In this case the functions A and B are given 

by 

It can then be seen that as i? ^ 0, vl and B approach the constant values 

where i_ is the lowest power in f{R). Similarly, as R ^ itcxD we have 

A^-i— and B ^'^-'^ 



2(i+ - 1) H ' 

where z+ is the highest power in f{R). In such limits these theories behave as if f{R) = R^, with the 
lowest power in (I84p dominating as i? ^ 0, and the highest power dominating as -R — > ±00. Even in 
the absence of a highest power in f{R) (that is, for an infinite power series) A and B can still approach 
constant values, a necessary condition to allow the existence of the vacuum dominated solutions, (jlOp . 
and matter dominated solutions, ()12|) . 
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The behaviour of R at late times, and hence the possible asymptotic behaviour of a(t), can now 
be determined from the behaviour of R in that limit. If the vacuum dominated evolution (jlOp is 
approached, as often is the case, then this behaviour can be read off from condition ()40p . For 1/2 < 
-B < 1 we have at late times that t ^ to, which corresponds to i? — > ±00 and, therefore, the 
asymptotic behaviour of a{t) being dominated by the highest power of R in the Lagrangian. This 
range of B corresponds to a power i?" with n > 2. Conversely, for i? > 1 or < — 1 we have that 
t ^ 00 at late times, so that -R — > 0, and there can exist an asymptote in which a{t) is dominated by 
the lowest power of R. These ranges of B correspond to R^ with n < 1/2. If a theory contains a low 
power of R in its Lagrangian there can then exist a vacuum dominated asymptote where i? — > 0, and 
in which the lowest power of R dominates. However, if that same theory contains high powers of R 
in its Lagrangian there can also exist an asymptote where R ^ 00, and in which the highest power 
of R dominates. Theories with multiple powers of R can therefore have their late-time asymptotics 
dominated by either the lowest or highest power of R. 

To illustrate these points further will now consider some example theories. 

5.1 Adding Higher Powers to General Relativity 

Consider adding to the Einstein-Hilbert action higher powers of R, such as R^ and R^. Such modifi- 
cations have often been considered as possible UV 'corrections' to gravity, motivated by attempts to 
construct a perturbatively re-normalisable quantum field theory of gravity. It is frequently assumed 
that these higher-order 'corrections' to gravity should dominate the earliest stages of the Universe's 
evolution. Let us investigate this possibility in the frame-work that has been constructed above. 
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Figure 7: The functions A{z) and B{z) for the theory f{R) = R + lO^R'^ + O.IR^. The dashed hue 
corresponds to A, and the solid line to B. 

Figure [7] shows the functions A[z) and B{z) when f{R) = R + lO^R^ + O.li?'^. The dashed line in 
this plot is A, and the solid line is B. One can see the linear term dominates at low R, the quadratic 
term at intermediate R, and the cubic term at high R. 

Consider a trajectory in the 'high-curvature' regime, z > 12, so that B ~ 2/3. In this case the 
stability of the various critical points, and their locations, are shown in Figure [21 The only expanding 
attractor is the vacuum dominated solution, (jlOp . which in this case corresponds to a ~ (t — ig)"^''- 
As this attractor is approached we have z' > 0, so trajectories are pushed further into the the high 
R regime, and toward an eventual big rip. For a trajectory at high-curvature to find its way down 
into the intermediate curvature regime, where R^ dominates, it must therefore be either collapsing, or 
possibly in the early transient stages before an asymptotic attractor is approached. 



23 



Now consider a trajectory in the B? dominated regime, —3 ^ z < h, where B ~ 1/2. The (x,y) 
plane looks very similar to Figure [2] in this case, but with the vacuum dominated points 5 and 6 
displaced a little. Again, the only expanding attractor is the vacuum dominated solution, which in 
this case corresponds to exponential expansion, (j26p . This point is then semi-stable, in the sense 
described in Appendix |Al Here we have AB^ > 0, so the vacuum dominated exponential expansion 
will eventually end with the point becoming a saddle in the phase space. Trajectories can then move 
to the higher or lower R regimes. 

The picture of a universe starting off expanding with the highest power in f{R) dominating, and 
moving successively down through lower powers of R until f r^ R is reached does not appear to be a 
generic situation. In fact, if a higher power of R (other than R^) dominates the gravitational dynamics, 
then it appears that expanding trajectories are generically forced to higher R, rather than lower. The 
R^ case is exceptional, and in this case expanding vacuum dominated trajectories do appear able to 
move to the lower R regime, although even in this case it does not seem generic. 



5.2 General Relativity as a Higher Power 

Consider now a gravitational Lagrangian with an Einstein-Hilbert term, and both higher and lower 
powers of R. Lower powers of R have been of interest recently, as they can correspond to cosmologies 
that accelerate at late times. 
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Figure 8: The functions A{z) and B{z) for the theory f{R) 
corresponds to A, and the solid line to B. 
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In Figure El we show the form of A{z) and B{z) for the theory f{R) = R+ W'^/R + W^R^. At 
high R the R^ term dominates, at intermediate R the Einstein-Hilbert term dominates, and at low 
R the inverse power, 1/R, dominates. As before, when R^ dominates the only expanding attractor 
at late times corresponds to increasing R. So again, only contracting solutions, or trajectories in the 
early transient stages, can make it down to the lower curvature regimes. A solution in the intermediate 
curvature regime, —3 ^ z ^ 3, is dominated here by the Einstein-Hilbert term, with S ~ 0. In this 
case the dynamical equations reduce to the usual Friedmann ones, and so R decreases for expanding 
solutions, allowing for a transition to the 1/R regime, z < —12, where i? ~ 2. 

The stability properties and location of critical points in the 1/R regime are shown in Figure [3l In 
this situation we again have that the only expanding attractor is the vacuum dominated solution, (llOp . 
In this case, however, A < 0, so this attractor corresponds to decreasing R. Trajectories asymptoting 
toward this point are then pushed further and further into the 1/R dominated regime. If we were 
to have included an even lower power of R, say l/R?, then expanding solutions heading toward the 
attractor solution would be pushed toward the regime in which this power dominates. 
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Expanding solutions in regimes dominated by the Einstein-Hilbert term, or lower powers of R, 
generically appear to move to lower and lower R as they continue to expand. The term R^ then acts 
as a watershed: expanding universes dominated by higher powers of R appear to evolve toward higher 
R, while expanding universes dominated by lower powers move to lower R. For the case of high R 
this corresponds to divergent expansion toward a big rip, and for low R it corresponds to (possibly 
accelerating) eternal expansion. 

5.3 Infinite Power Series 

Having considered the case of single higher or lower powers dominating the gravitational Lagrangian, 
let us now consider infinite power series. There are, of course, very many functional forms that one 
may choose for f{R) to illustrate the case of an infinite power series. Here we consider f{R) = 
{R-\- 100i?^)e^'^ . The corresponding functions A[z) and B{z) are shown in Figure[9l It can be seen 
that there are three different regimes: low i? at z < —5, where the Einstein-Hilbert term dominates 
with i? ~ 0, intermediate i? at < z < 10, where R"^ dominates with B ~ 2/3, and high R at z > 15, 
where the exponential dominates with B ~ 1. 

A(z),B(z) 
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Figure 9: The functions A{z) and B{z) for the theory f{R) = {R + 100R^)e^^~^^. The dashed hue 
corresponds to A, and the solid line to B. 

The high R regime now corresponds to 5 ~ 1, which is one of the special cases in Appendix lAl 
The vacuum dominated point 5 is again the only expanding attractor, and now this point corresponds 
to exponential growth of the scale factor, as in ()26p . This point is then stable, as Az ^ 0, but with 
z' = 0. While it can be the case that single higher powers dominate at intermediate regimes (i?^ in 
the region < z < 10, here), the solutions heading toward the expanding attractor in such regimes are 
forced to higher R. The highest R regime here, however, is at B ~ 1. The late-time attractor then no 
longer corresponds to a big rip, but to exponential expansion, with R ^constant, instead of oo. The 
behaviour of infinite power series in R is therefore qualitatively different to the finite series considered 
previously. As before, only collapsing solutions, or possibly those in the early transient stages of their 
evolution, appear able to make it down to the low curvature regime. 

6 Discussion 

We have considered in this paper the evolution of spatially flat FRW universes governed by f{R) 
theories of gravity. The Friedmann equations ([6])- ([8]) were transformed into an autonomous system of 
first-order differential equations, and a dynamical systems analysis was performed. The location and 
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stability of all critical points in the phase space were found, for both vacuum and perfect fluid cos- 
mologies, and for general f{R)- It was shown that the simple power-law solutions, given by equations 
(|10 p -(jl3 p . often act as the early and late-time asymptotes of the general solution. 

The general behaviour of f{R) FRW cosmologies is complicated. The phase space of solutions is 
often divided into sub-spaces by invariant manifolds, through which the trajectories describing the 
general solution cannot pass. The asymptotic past of general solutions can contain points of inflexion, 
or big-bang singularities that can be approached in different ways. Similarly, future behaviour can be 
seen to be able to asymptote toward matter dominated expansion, (J12p . vacuum domination, (jlOp . or 
various other forms. Whatsmore, the three dimensional phase space of solutions, in the presence of a 
perfect fluid, allows for the possible existence of strange attractors, and chaotic behaviour. 

Nevertheless, despite the complicated behaviour exhibited by the general solutions, it is still pos- 
sible to make statements about the effects of modifying the Einstein-Hilbert action to more general 
functions of R. It can be said that theories that contain lower powers of R often have a stable asymp- 
tote that corresponds to the expanding vacuum dominated solution, (jlOp . of that lowest power, and 
that therefore behave as if governed by a gravitational Lagrangian of the form ([2]) at late times. How- 
ever, if a theory contains any higher powers of R then there is also often a stable asymptote that 
corresponds to the expanding vacuum dominated solution, (llOp . of that highest power, and that there- 
fore behaves as if governed by a Lagrangian of the form ([3]) at late times. Theories containing both 
higher and lower powers of R can then asymptote, at late times, to regimes in which either the lowest 
or highest powers or R dominate. In either case the consequent evolution is that of a gravitational 
Lagrangian dominated by a single power of R. 

Which solutions asymptote to high R domination, and which to low R domination, depends on the 
initial conditions, and the form of f{R)- Using illustrative examples we have shown that if a power 
of R higher than R^ dominates at some point then the generic behaviour of expanding solutions is 
to higher R. If a power of R lower than R^ dominates, then the trend is to lower R. A term R^ 
in the Lagrangian is then a special case; if it dominates then the expanding attractor corresponds 
to exponential growth, and acts as a separatrix between the higher or lower powers of R dominating 
the future dynamics of the Universe. Those solutions with low R dominating at late-times expand 
eternally, while those with high R dominating approach either a big-rip singularity, or exponential 
expansion. Big-rips often occur if there exists a single power of R that dominates at late times, 
and exponential expansion occurs if i? — > 1, as is the case for some infinite power series, such as 
/ ~ exp{i?}. This late-time acceleration does not appear to be a good candidate for the apparent 
acceleration we observe, however, as it cannot follow from a period of Einstein-Hilbert domination in 
which i? — > 0. 

The asymptotic past of the general solutions is similarly complicated. A variety of behaviours 
seems possible, including big-bang singularities and bounces, where the scale factor reaches a non-zero 
minimum. Big-bang singularities are often approached with the scale-factor behaving as in the Tolman 
solution, (J13p . It is interesting to note that even trajectories which undergo an early period of inflation 
(such as those solutions in which a power of R^ dominates at some point) do not generically follow 
such expansion indefinitely into the past, but rather have a big bang or bounce at some point in their 
past, prior to the onset of infiation. 

While the behaviour found here is quite complicated, with numerous different asymptotes possible, 
it is certainly not the most general case that one may consider. We have limited ourselves here 
to spatially fiat FRW universes. Relinquishing the criterion of spatial flatness would lead to more 
complicated behaviours still, and one may also consider inhomogeneous and/or anisotropic cosmologies, 
or cosmologies with multiple fluids. It is not clear whether or not the behaviour identifled above would 
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hold in these more general cases or not. What does seem clear, however, is that the simple picture 
of the highest powers of R dominating at early times, and lower powers dominating at late times, is 
unlikely to be accurate. 

A Stability of 'Vacuum Dominated' Solutions when B = 1/2 or 1 

When B = 1/2 or 1 the vacuum dominated solutions take an exponential form a ~ e'^'*"*''^ instead of 
the usual power-law form, (llOh . In this case the parameter z = In |i?| does not diverge to ±00, and so 
the stability analysis must be modified from the power-law case. We will investigate the stability of 
these solutions here. 

A.l Vacuum Cosmologies 

In the vacuum cosmologies, considered in section 3, points 5 and 6 correspond to exponential expansion 
when B — > 1/2 or 1. Consider first the case B — > 1/2. Perturbing 9 ^> 9 + 56 and z ^ z + 6z we then 
have that the evolution equations ()19p and (120p are, to linear order. 



S0' = ^^69 ± ^B^Jz (89) 

2 \/3 

Sz' = ±aJ^69, (90) 

V o 

where the upper branch corresponds to point 5, and the lower branch to point 6. It can now be seen 
that the eigenvalues to the equations 5z' = Xi6z and 59' = Xi69 are given by the two roots of 

.2 , V3Q, 4 



Xf ± ^Xi - -AB^.Q = 0. (91) 

On recognising that these solutions have Q = 1^ \i can be seen that if AB^z > then both points are 
saddles. Alternatively, if AB^ < then point 5 is stable, while point 6 is unstable. We will also be 
interested in the case when B^z — 0. Point 5 is then stable in 9, while Sz' = 0. Small fluctuations (as 
one would expect to occur in a thermal de Sitter space such as this) will then gradually shift the value 
of z. Unless the theory has B^ = for all z, there will then come a time at which B^ ^ 0, after which 
the point becomes an attractor proper, or a saddle, depending on the sign of AB^- We will refer to 
this situation of temporary stability, followed by saddle behaviour, as 'semi-stable'. 

Now consider the case B ^ 1. In this case it is instructive to use the definitions of A and B to see 
that 

For Bz 7^ —1 we then have that B ^ 1 corresponds to A — > 0. Perturbing 9 and z, as before, we then 
find that the linearised evolution equations become 




(93) 
(94) 
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so that the relevant eigenvalues are given as the roots of 

A? ± y|(2 - A,)A, - ^A^Q = 0. (95) 

Again, Q = 1, so the condition for points 5 and 6 to be saddles is now A^^ > 0. If A^z < then point 
5 is an attractor, and point 6 a repellor. Now, if A^z = 0, then point 5 can be semi-stable, depending 
on the sign of Az ■ 

A. 2 Perfect Fluid Cosmologies 

In the perfect fluid cosmologies we also find that the vacuum dominated solution corresponds to 
exponential expansion when B -^ 1/2 or 1. Consider first the case B — > 1/2. In this case a perturbative 
expansion in </> and z gives the evolution equations (j50p . (j58p and (j59p . to linear order, as 

(96) 

5z' = ±yjp^ (97) 

(1-r) 

/ = ±^, (98) 

where the upper sign is for point 5 and the lower sign for point 6. We are now looking for eigenvalues, 
Aj, such that 

69' = -^^69 and 6z' = -^^6z. (99) 

(1 — r) (1 — r) 

If both of these eigenvalues are negative, and r' > 0, then the point is stable. If they are both positive, 
and r' < 0, then the point is unstable. Any other combination we will call a saddle point. The values 
of Aj can be seen to be given by the roots of 

Xj ± 3Ai - 16AB^z = 0. (100) 

It can then be seen that for AB z > both points are saddles, while for AB z < point 5 is stable 
and point 6 is unstable. ABz = again corresponds to point 5 being semi-stable. The stability 
properties of points 5 and 6 in the presence of a perfect fluid when B — > 1/2 are therefore identical to 
the corresponding points in the vacuum cosmologies, considered above. 

Finally, consider the case i? — > 1. Equation (j92p shows that in this limit we again have ^ ^ 0. 
The perturbed evolution equations are then given as 

^^^^jw^z-m (,,,) 

V2{l-r) 

^ , \/2Az6z 

6z' = ± ; '\ 102 

(1-r) 

r' = ±^(l + 37), (103) 

where upper signs correspond to point 5, and lower signs to point 6. The eigenvalues, Aj, are now 
given by the roots of 

Af ± V2(2 - A^z)\i - 4A^z = 0. (104) 

28 



For A^z > both of points 5 and 6 are saddles, while for A^z < point 5 is stable and point 6 
is unstable. A^z = gives point 5 as semi-stable. Again, the stability properties of these points is 
identical to the analogous points in the vacuum cosmologies. 

Acknowledgements 

I would like to thank John Barrow for suggestions, and to acknowledge the support of Jesus College, 
Oxford. 



References 



[9 

[lo; 

[11 

[12 

[13; 

[14 
[15 

[16; 

[ir 
[is: 

[19 

[2o; 

[21 



H. Buchdahl, J. Phys. A 12, 1229 (1979). 

R. Kerner, Gen. Rel. Grav. 14, 453 (1982). 

J. D. Barrow and A. C. Ottewill, J. Phys. A 16, 2757 (1983). 

G. Magnano, M. Ferraris and M. Francaviglia, Gen. Rel. Grav. 19, 465 (1987). 

E. Pechlaner and R. Sexl, Comm. Math. Phys. 2, 165 (1966). 

M. Gasperini, M. Maggiore and G. Veneziano, Nucl. Phys. B 494, 315 (1997). 

K. A. Meissner, Phys. Lett. B 392, 298 (1997). 

S. M. Carroll, A. De Felice, V. Duvvuri, D. A. Easson, M. Trodden, and M. S. Turner, Phys. 
Rev. D 71, 063513 (2005). 

S. Nojiri and S. D. Odintsov, Phys. Rev. D 68, 123512 (2003). 

A. Berkin, Phys. Rev. D 44, 1020 (1991). 

E. Bruning, D. Coule and C. Xu, Gen. Rel. Grav. 26, 1197 (1994). 

J. D. Barrow and S. Hervik, Phys. Rev. D 73, 023007 (2006). 

J. D. Barrow and J. Middleton, Phys. Rev. D 75, 123515 (2007). 

J. D. Barrow and T. Clifton, Class. Quant. Grav. 23, LI (2006). 

T. Clifton and J. D. Barrow, Class. Quant. Grav. 23, 2951 (2006). 



T. P. Sotiriou and V. Faraoni, arXiv:0805.1726 (2008). 

T. Clifton and J. D. Barrow, Phys. Rev D 72, 103005 (2005). 

S. Carloni, P. K. S. Dunsby, S. Capoziello and A. Troisi, Class. Quant. Grav. 22, 4839 (2005). 

M. Abdelwahab, S. Carloni and P. K. S. Dunsby, Class. Quant. Grav. 25, 135002 (2007). 

J. A. Leach, S. Carloni and P. K. S. Dunsby, Class. Quant. Grav. 23, 4915 (2006). 

N. Goheer, J. A. Leach and P. K. S. Dunsby, Class. Quant. Grav. 24, 5689 (2007). 



29 



[22] J. D. Barrow and S. Hervik, Phys. Rev. D 74, 124017 (2006). 

[23] T. Clifton, Class. Quant. Grav. 24, 5073 (2007). 

[24] L. Amendola, R. Gannouji, D. Polarski and S. Tsujikawa, Phys. Rev. D 75, 083504 (2007). 



[25] S. Carloni, A. Troisi and P. K. S. Dunsby, arXiv:0706. 0"452 (2007). 
[26] T. Clifton, Class. Quant. Grav. 23, 7445 (2006). 
[27] T. Clifton, Phys. Rev. D 77, 024041 (2008). 



30 



